% Copyright (C) 2014-2022 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer
%
% This is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This code is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% For a copy of the licencse,
% see <http://www.gnu.org/licenses/>.

% Create Figure 8: Balanced budget government spending shock
clear all
% close all
addpath('../Auxiliary_Files')
if ~isfolder('results')
    mkdir('.','results');
end
if ~isfolder('Figures')
    mkdir('.','Figures');
end
if ~isfolder('LP_coefficients')
    mkdir('.','LP_coefficients');
end


for sign_iter=1:2
    for shock_iter=1:2
        store_coefficients_T_G(sign_iter,shock_iter)
    end
end

rng(1);
ptileVEC = [0.16 0.50 0.84]; % Percentiles for posterior distributions
% ptileVEC = [0.05 0.50 0.95]; % Percentiles for posterior distributions
transp90 = 0.25;
fontsize = 6;
linW = 2;

g.hike=load('LP_coefficients\ghike.mat');
g.cut=load('LP_coefficients\gcut.mat');
t.hike=load('LP_coefficients\thike.mat');
t.cut=load('LP_coefficients\tcut.mat');

Horizon=size(g.cut.theta_mat,1);
n_vars=size(g.cut.theta_mat,2);
n_restriction_periods=4;

shock_names={'cut','hike',};
var_names={'G','Y','FX (down is apprecation)','T'};
FX_pos=3;

n_draws=10000;
hh=figure('Name','Balanced Budget Shock');
for shock_iter=1:2
    IRF_mat.(shock_names{shock_iter})=zeros([n_vars,Horizon,n_draws]);
    for draw_iter=1:n_draws
        g.(shock_names{shock_iter}).theta_mat_draw=g.(shock_names{shock_iter}).theta_mat+g.(shock_names{shock_iter}).se_mat.*randn(size(g.(shock_names{shock_iter}).theta_mat));
        t.(shock_names{shock_iter}).theta_mat_draw=t.(shock_names{shock_iter}).theta_mat+t.(shock_names{shock_iter}).se_mat.*randn(size(t.(shock_names{shock_iter}).theta_mat));
        if strcmp(shock_names{shock_iter},'hike')
            target=[ones(1,n_restriction_periods); ones(1,n_restriction_periods)]; % raise revenue by relative ratios to keep budget balanced, G by 1
        elseif strcmp(shock_names{shock_iter},'cut')
            target=-[ones(1,n_restriction_periods); ones(1,n_restriction_periods)]; % raise revenue by relative ratios to keep budget balanced, G by 1
        else
            error('Not defined')
        end
        shock_mat.(shock_names{shock_iter})=zeros(2,9);

        %get shock matrices
        A=[g.(shock_names{shock_iter}).theta_mat_draw(1,1) t.(shock_names{shock_iter}).theta_mat_draw(1,1);
            g.(shock_names{shock_iter}).theta_mat_draw(1,end) t.(shock_names{shock_iter}).theta_mat_draw(1,end)];
        for h_iter=0:n_restriction_periods-1
            for lag_iter=1+(h_iter-1):-1:1
                target(:,1+h_iter)=target(:,1+h_iter)-[g.(shock_names{shock_iter}).theta_mat_draw(1+lag_iter,1)   t.(shock_names{shock_iter}).theta_mat_draw(1+lag_iter,1)
                                                       g.(shock_names{shock_iter}).theta_mat_draw(1+lag_iter,end) t.(shock_names{shock_iter}).theta_mat_draw(1+lag_iter,end)]*shock_mat.(shock_names{shock_iter})(:,1+h_iter-lag_iter);
            end
            shock_mat.(shock_names{shock_iter})(:,1+h_iter)=A\target(:,1+h_iter);
        end


        %compute IRFs
        for h_iter=0:Horizon-1
            for lag_iter=0:h_iter
                IRF_mat.(shock_names{shock_iter})(:,1+h_iter,draw_iter)=IRF_mat.(shock_names{shock_iter})(:,1+h_iter,draw_iter)+[g.(shock_names{shock_iter}).theta_mat_draw(1+lag_iter,:)' t.(shock_names{shock_iter}).theta_mat_draw(1+lag_iter,:)']*shock_mat.(shock_names{shock_iter})(:,1+h_iter-lag_iter);
            end
        end
    end
    temp=quantile(IRF_mat.(shock_names{shock_iter}),ptileVEC,3);
    temp(FX_pos,:,:)=-temp(FX_pos,:,:); %flip FX due to notation
    for var_iter=1:length(var_names)-1% do not do taxes
        subplot(3,2,shock_iter+(var_iter-1)*2)
        hold on
        a = squeeze(temp(var_iter,:,1));
        b = squeeze(temp(var_iter,:,3));
        ha1_90 = area(0:Horizon-1,[a; b-a]','FaceColor',[153/255 204/255 1],'EdgeColor','none','ShowBaseLine','off');
        set(ha1_90(1), 'FaceColor', 'none') % this makes the bottom area invisible
        set(ha1_90, 'LineStyle', '-')
        hold on
        aa=plot(0:1:Horizon-1,squeeze(temp(var_iter,:,2)),'b','LineWidth',1);
        if strcmp(var_names(var_iter),'G')
            tax_pos=strmatch('T',var_names,'exact');
            a = squeeze(temp(tax_pos,:,1));
            b = squeeze(temp(tax_pos,:,3));
            ha1_90 = area(0:Horizon-1,[a; b-a]','FaceColor',[1 0.6 0.6],'EdgeColor','none','ShowBaseLine','off');
            set(ha1_90(1), 'FaceColor', 'none') % this makes the bottom area invisible
            set(ha1_90, 'LineStyle', '-')
            ha1_90(2).FaceAlpha = 0.5;
            hold on
            bb=plot(0:1:Horizon-1,squeeze(temp(tax_pos,:,2)),'r--','LineWidth',1);
            if shock_iter==1
                ll=legend([aa,bb],'G','T');
                ll.Location='NorthWest';
                ll.Box='off';
            end
        end
        hline(0,'k:');
        %         set(gca,'YTick',-0.5:0.5:2.5)
        set(gca,'XTick',0:4:Horizon-1)
        xlim([0 Horizon-1])
        box on;set(gca,'xTick',0:Horizon-1,'Layer','top','FontSize',fontsize);
        if var_iter==1
            title([shock_names{shock_iter} ': ' var_names{var_iter}])
        else
            title(var_names{var_iter})
        end
    end
end

print(['Figures/LP_balanced_budget_', num2str(n_restriction_periods)],'-dpdf')
print(['Figures/LP_balanced_budget_', num2str(n_restriction_periods)],'-depsc2')
saveas(hh,['./Figures/LP_balanced_budget'])